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Abstract 

We apply a probabilistic approach to study the computational complexity of ana- 
log computers which solve linear programming problems. We analyze numerically 
various ensembles of linear programming problems and obtain, for each of these en- 
sembles, the probability distribution functions of certain quantities which measure 
the computational complexity, known as the convergence rate, the barrier and the 
computation time. We find that in the limit of very large problems these proba- 
bility distributions are universal scaling functions. In other words, the probability 
distribution function for each of these three quantities becomes, in the limit of large 
problem size, a function of a single scaling variable, which is a certain composition 
of the quantity in question and the size of the system. Moreover, various ensembles 
studied seem to lead essentially to the same scaling functions, which depend only on 
the variance of the ensemble. These results extend analytical and numerical results 
obtained recently for the Gaussian ensemble, and support the conjecture that these 
scaling functions are universal. 

PACS numbers: 5.45-a, 89.79-hc, 89.75.D 



1 Introduction 

Digital computers are part of our present civilization. There are, however, other devices 
that are capable of computation. These are analog computers, that are ubiquitious com- 
putational tools. The most relevant examples of analog computers are VLSI devices im- 
plementing neural networks or neuromorphic systems p'^, whose structure is directly 
motivated by the workings of the brain. Various processes taking place in living cells can 
be considered as analog computation |H] as well. 

An analog computer is essentially a physical device that performs computation, evolving 
in continuous time and phase space. It is useful to model its evolution in phase space by 
dynamical systems (DS) |5, the way classical systems such as particles moving in a potential 
(or electric circuits), are modeled. For example, there are dynamical systems (described 
by ordinary differential equations) that are used to solve computational problems [HI El EI ■ 
This description makes a large set of analytical tools and physical intuition, developed for 
dynamical systems, applicable to the analysis of analog computers. 

In contrast, the evolution of a digital computer is discrete both in its phase space and 
in time. Consequently, the standard theory of computation and computational complexity 
[SI deals with computation in discrete time and phase space, and is inadequate for the 
description of analog computers. The analysis of computation by analog devices requiers a 
theory that is valid in continuous time and phase space. 
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Since the systems in question are physical systems, the computation time is the time 
required for a system to reach the vicinity of an attractor (a stable fixed point in the present 
work) combined with the time required to verify that it indeed reached this vicinity. This 
time is the elapsed time measured by a clock, contrary to standard computation theory, 
where it is the number of discrete steps. 

In the exploration of physical systems, it is sometimes much easier to study statistical 
ensembles of systems, estimating their typical behavior using statistical methods [Hl llOlITT] . 
In jl2l I13j a statistical theory was used to calculate the computational complexity of a 
standard representative problem, namely Linear Programming (LP), as solved by a DS. 

A framework for computing with DS that converge exponentially to fixed points was 
proposed in |14j . For such systems it is natural to consider the attracting fixed point as 
the output. The input can be modeled in various ways. One possible choice is the initial 
condition. This is appropriate when the aim of the computation is to decide to which 
attractor out of many possible ones the system flows 15 . Here, in |121 113j . as well as in 
[14j . the parameters on which the DS depends (e.g., the parameters appearing in the vector 
field F in (P)) are the input. 

The basic entity of the computational model is a dynamical system 0], that may be 
defined by a set of Ordinary Differential Equations (ODEs) 



where x is an n-dimensional vector, and F is an n-dimensional smooth vector field, which 
converges exponentially to a fixed point. Eq. solves a computational problem as follows: 
Given an instance of the problem, the parameters of the vector field F are set (i.e., the 
input is read), and it is started from some pre-determined initial condition. The result of 
the computation is then deduced from the fixed point that the system approaches. 

In our model we assume we have a physical implementation of the flow equation (^. 
Thus, the vector field F need not be computed, and the computation time is determined by 
the convergence time to the attractive fixed point. In other words, the time of flow to the 
vicinity of the attractor is a good measure of complexity, namely the computational effort, 
for the class of continuous dynamical systems introduced above jl4j . 

In this paper, as in |12| I13j . we will study a specific algorithm for the solution of the 
LP problem ^H]- We will consider real-continuous inputs, as the ones found in physical 
experiments, and that are studied in the BSS model J7j) ^ '^^11 as integer- valued inputs. 
For computational models defined on the real numbers, worst case behavior, that is tra- 
ditionally studied in computer science, can be ill-defined and lead to infinite computation 
times, in particular, for some methods for solving LP |17| I18j. Therefore, we compute the 
distribution of computation times for a probabilistic model of LP instances with various 
distributions of the data like in |19[ I2nj . Ill-defined instances constitute a set of measure 
zero in our countinuous probability ensembles, and need not be concerned about. In the 
discrete probability ensembles, we treat them by appropriate regularization. 

The computational complexity of the method presented in E| and discussed here 
is 0{nlogn), compared to 0{n^'^ logn) found for standard interior point methods |21| . 
The basic reason is that for standard methods (such as interior point methods), the major 
component of the complexity of each iteration is 0{n'^) due to matrix decomposition and 
inversion of the constraint matrix, while here, because of its analog nature, the system just 
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flows according to its equations of motion (which need not be computed). 

Since we consider the evolution of a vector field, our model is inherently parallel. There- 
fore, to make the analog vs. digital comparison entirely fair, we should compare the com- 
plexity of our method to that of the best parallel algorithm. The latter can reduce the 
0{n^) time needed for matrix decomposition/inversion to poly logarithmic time (for well- 
posed problems), at the cost of ©(n^-^) processors [211, while our system of equations 
uses only 0{n) variables. 

The main result of jl2l I13j . in which LP problems were drawn from the Gaussian 
distribution of the parameters of F (namely, the constraints and cost function in ((21)), was 
that the distribution functions of various quantities that characterize the computational 
complexity, were found to be scaling functions in the limit of LP problems of large size. In 
particular, it was found that these distribution functions depend on the various parameters 
only via specific combinations, namely, the scaling variables. Such behavior is analogous 
to the situation found for the central limit theorem, for critical phenomena |22] and for 
Anderson localization [2^, in spite of the very different nature of these problems. It was 
demonstrated in ^r71ll3j how for the implementation of the LP problem on a physical device, 
methods used in theoretical physics enable to describe the distribution of computation times 
in a simple and physically transparent form. Based on experience with certain universality 
properties of rectangular and chiral random matrix models j2Sl) it was conjectured in |121 
E] that some universality for computational problems should be expected and should be 
explored. That is, the scaling properties that were found for the Gaussian distributions 
should hold also for other distributions. In particular, some specific questions were raised 
in jl2l I13j : Is the Gaussian nature of the ensemble unimportant in analogy with UHl? Are 
there universality classes [2Sj of analog computational problems, and if they exist, what are 
they? 

Thus, we extend the earlier analysis |12[ ll3j of the Gaussian distribution to other prob- 
ability distributions of LP problems, and demonstrate numerically that the distribution 
functions of various quantities that characterize the computational complexity of the ana- 
log computer which solves LP problems are indeed universal scaling functions, in the limit 
of large systems. These universal functions depend upon the original probability ensemble 
of inputs only via the scaling variables, that are proportional to the ones found for the 
Gaussian distribution. For some distributions of LP problems the scaling variables are even 
identical (not just proportional) to the ones found for the Gaussian distribution. For other 
distributions, on the other hand, where some of the parameters defining an LP problem 
may vanish at random (the so-called diluted ensembles in Section 4) , either the convergence 
to universality is much slower, or universality is only approximate. 

The distribution of constraints and cost function of the LP problems that are used 
in practice is not known. Therefore, the universality of the distribution functions of the 
computation time and other quantities related to computational complexity is of great 
importance. It would imply that it may hold also for the distributions of the LP problems 
solved in applications. In this paper we demonstrate numerically that for several probability 
distributions universality is satisfied, providing support for the conjecture that it holds in 
general. 

This paper is organized as follows: In the next section we briefly review the dynamical 
system which solves LP problems. In section 3 we summarize the scaling results of jl2l I13j 
for the Gaussian ensemble. In section 4 we present our numerical results for the distribution 
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functions of various quantities that characterize the computational complexity of the analog 
computer for non-Gaussian probability ensembles. In section 5 we demonstrate that these 
distributions are indeed universal scaling functions. Finally, in section 6 we discuss the 
significance of our results and also pose some open problems. 
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2 A dynamical flow for linear programming 



Linear programming is a P-complete problem (8j|, i.e. it is representative of all problems 
that can be solved in polynomial time. The standard form of LP is to find 

max{c^x : x e , Ax = b,x > 0} (2) 

where c G M", 6 G M™", A G R™^" and m < n. The set generated by the constraints in Q 
is a polyheder. If a bounded optimal solution exists, it is obtained at one of its vertices. 
The vector defining this optimal vertex can be decomposed (in an appropriate basis) in the 
form X = {xj\f,XB) where xj\f = is an n — m component vector, while = B~^b > is 
an m component vector, and B is the m x m matrix whose columns are the columns of A 
with indices identical to the ones of xg. Similarly, we decompose A = (N,B). 

A flow of the form converging to the optimal vertex, introduced by Faybusovich |Sj 
will be studied here. Its vector fleld F is a projection of the gradient of the cost function 
c^x onto the constraint set, relative to a Riemannian metric which enforces the positivity 
constraints x > [0]. It is given by 

F{x) = [X - XA^{AXA'^)-^AX] c , (3) 

where X is the diagonal matrix Diag(xi . . . The nm + n entries of A and c, namely, the 
parameters of the vector field F, constitute the input; as in other models of computation, we 
ignore the time it takes to "load" the input, since this step does not refiect the complexity of 
the computation being performed, either in analog or digital computation. It was shown in 
jl6j that the flow equations given by and ^ are, in fact, part of a system of Hamiltonian 
equations of motion of a completely integrable system of a Toda type. Therefore, like the 
Toda system, it is integrable with the formal solution jB] 

x,{t) = xM exp I -A,t + V a,, log ^2±nz!!igl 1 (4) 



/ . "J 



m(0) 



(i = 1, . . . ,n — m), that describes the time evolution of the n — m independent variables 
xj\f{t), in terms of the variables xs{t). In Q Xi{0) and Xj+n-m{0) are components of the 
initial condition, Xjj^n-m{i) are the xs components of the solution, aji = —{B^^N)ji is an 
m X (n — m) matrix, while 

m 

Ai = -Ci-Y^ Cjaji ■ (5) 
i=i 

For the decomposition x = {xj\f, xq) used for the optimal vertex Aj>0 i = l,...,n — m, 
and Xf/(t) converges to 0, while converges to x* = B~^b. Note that the analytical 

solution is only a formal one, and does not provide an answer to the LP instance, since 
the Aj depend on the partition of A, and only relative to a partition corresponding to a 
maximum vertex are all the Aj positive. 

The second term in the exponent in Q), when it is positive, is a kind of "barrier": Ajt 
must be larger than the barrier before Xi can decrease to zero. In the following we ignore 
the contribution of the initial condition and denote the value of this term in the infinite 
time limit by 

m 

A = II"iilog^i+n-m- (6) 
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Note that although one (or more) of the may vanish, in the probabiUstic ensemble 

studied here, such an event is of measure zero for the continuous probability distributions, as 
well as for the regularized discretized ones (see (US))), and therefore should not be considered. 
In order for x{t) to be close to the maximum vertex we must have Xi{t) < e for i = 
1, . . . , n — m for some small positive e, namely exp(— Ajt + < e , for i = 1, . . . , n — m. 
Therefore we consider 

f Pi \loge\\ 

as the computation time. We denote 

Amin = minAj, /?max = max /Jj . (8) 

The Aj can be arbitrarily small when the inputs are real numbers, but in the probabilistic 
model, "bad" instances, resulting in computation taking arbitrarily long time, are rare as 
is clear from^ Q. 



^Strictly speaking, Q was derived in |12II13| for a probabilistic model in which the components of {A, b, c) 
were independent identically distributed Gaussian random variables. However, one of the main points of 
this paper is that Q is valid for a broad class of probability distributions of LP problems. 
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3 Probability distributions and scaling 



Consider an ensemble of LP problems in which the components of {A, b, c) are independent 
identically distributed (i.i.d.) random variables taken from various even distributions, with 
mean and bounded variance. In a probabilistic model of LP instances, Amim /3max and 
T are random variables. Since the expression for Aj, equation ©, is independent of b, its 
distribution is independent of b. For a given realization of A and c, with a partition of A 
into (A^, B) such that Aj > 0, there exists^ a vector b such that the resulting polyheder has a 
bounded optimal solution. Since b in our probabilistic model is independent of A we obtain: 
7^(Aniin < A|A„jin > 0, LP instance has a bounded maximum vertex) = 'P(Ainin < AjAi^in > 
0). 

In jl2l I13j . the components of A, b, and c were taken from the Gaussian distribution 
(see, e.g., Eqs. (12-18) in jl2j ) with zero mean and variance o"^, that was taken as unity in 
the numerical calculations. It was found analytically, in the large (n, m) limit, that the 
probability P(Amm < A|A™.„ > 0) = ( A) is of the scaling form 

^(".'^)(A) = l-e^ierfc(xA) = .^(xa). (9) 

with the scaling variable 

XA{n,m) = ^ 1 ^ . (10) 



vr \m 

The scaling function contains all asymptotic information on A. The probability density 
function derived from J^{x/\) is very wide and does not have a finite variance. Also the 
average of l/x/\ diverges. 

The amazing point is that in the limit of large m and n, the probability distribution 
of Amin depends on the variables m, n and A only via the scaling variable xa- For future 
reference it is convenient to write xa in the form 

XA = CL^^\n/ 111)^/711 A , (11) 

with 

a'i\n/m) = ^(^-l)l, (12) 



'TT \m J a 

where the superscript refers to the Gaussian distribution. If the limit of infinite m and 
n is taken, so that njm is fixed, af is constant. It was verified numerically that for the 
Gaussian ensemble © was a good approximation already for m = 20, and n = 40. 

The existence of scaling functions like © for the barrier f3max^ that is the maximum of 
the /3j defined by © and for T defined by (assuming that e is not too small so that the 
first term in Q dominates) was verified numerically for the Gaussian distribution |12| ll.S) . 
In particular for fixed m/n, we found that 

V{l/f3^ax < 1//3) = -^J^^LlV/?) ^ ^i/fsixp) (13) 



^The existence of 6 is guaranteed by the fact that the various probability distributions are even. See |12| . 
Lemma 3.1. The latter was proved in under the assumption that the components of {A,b,c) were i.i.d. 
Gaussian variables, but the proof extends trivially to the class of probability ensembles of the type specified 
above. 
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and 

V{l/T < l/t) = ^ T^,t{xt). (14) 

The corresponding scaling variables are 

xp = af{n/m)'^ (15) 

and 

(g) mlogm 
XT = a>f{n/m) . [lb) 

Since the distribution functions ()13() and 1)141) are not known analytically, and since n = 2m 
was taken in the numerical investigations in J2], we can set arbitrarily a^^(2) = ai^\2) = 1. 
The scaling functions @, and (|14() imply the asymptotic behavior 

l/Amm~\/w^, /3max~m, i ~ m log m (17) 
with "high probability" [T^ITH] . 
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4 Non-Gaussian distributions 



In this section we present the results of our numerical calculations of the distribution 
functions of Amm, l/zSmax, and 1/T for various probability distributions of A, b, and c. 
For this purpose we generated full LP instances {A, b, c) with the probability distribution 
in question. For each instance the LP problem was solved using the linear programming 
solver of MatLab. Only instances with a bounded optimal solution were kept, and Amin 
was computed relative to the optimal partition and optimality was verified by checking that 
Amin > 0. Using the sampled instances we obtain an estimate of jr{«,"^)(A) = P(A™„ < 
A\Amin > 0), and of the corresponding cumulative distribution functions of the barrier 
/?max and the computation time T. 

The solution of the LP problem is used here in order to identify the optimal partition 
of A into B and N. This enables one to compute Amin, Pmax, and t from ©, © and 
(jHJ, and the distributions 

Pa(A) = V{Amin < A|A„i„ > 0) , (18) 

Pf3il//3)=V{l/Pma.<l/P) (19) 

and 

PT{l/t)=V{l/T <l/t) (20) 

are obtained. 

It is convenient at first to keep n/m fixed^ and to compute these distributions as func- 
tions of the corresponding scaling variables 

x'a = (21) 

x'f3 = m/P (22) 

and 

x'rp = mlogm/t (23) 

These are proportional to XA,a;/3 and xt defined by (fTU]) . ((T3)) and H16() . We turn now to 
explore in some detail, various distributions. 



4.1 The bimodal distribution 

In this case, the various elements of A, b and c take only the values +1 or —1 with probability 
1/2 each, namely, 

P{y) = \[6{y-l)+6{y + l)]. (24) 

The mean of this distribution is 0, and its variance is 1. One problem associated with the 
discrete ensemble ()24p is the finite probability to draw a degenerate LP problem. (Note 
that such degenerate, ill-defined LP problems, comprize a set of zero measure in continuous 
ensembles such as the Gaussian ensemble.) In order to avoid these degenerate solutions, 

^In fact, in all our numerical simulations we kept n/m — 2 fixed. 
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Figure 1: Va is plotted as a function of A, for the bimodal distribution (for n = 2m). 
The number of instances used in the simulation is 39732 for the m = 20 case, 46583 for the 
m = 30 case and 47169 for the m = 40 case. The number of converging instances for each 
case is 5000. 

we introduce a "regularization" by which each matrix element Aij, chosen at random from 
the ensemble (|2^ . is multiplied by 



where e is a small regularization parameter. (The other entries 6j and Cj take the values 
chosen at random with the probability (|24() without any regularization.) In the numerical 
calculations we take e = 10~^. (Note that the regularization (|25j) slightly splits the identical 
unregulated probability distributions of the independent matrix elements Aij. For small 
e we expect this fact to have a negligible effect on the scaling behavior.) In Figs. ^-El 
the distribution functions ((TH|) , ((T^ , and (|^n|) are plotted first as functions of the unsealed 
variables A, 1//3, and l/t, followed by the corresponding plots as functions of the scaled 
variables ()21() . (|22() . and ()23() . respectively, for various values of m while n = 2m. Note that 
in terms of the scaled variables (PT|) -(P ^ . the data is found to collapse to one distribution. 



/. . = 1 + [i + 2(i - l)]e 



(25) 
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Va is plotted as a function of x'^ for the bimodal distribution for the instances 
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Figure 3: Vfs is plotted as a function of for the instances of Fig. ^ 
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Figure 4: is plotted as a function of x'^ for the instances of Fig. ^ 
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Figure 5: Vt is plotted as a function of 1/t for the instances of Fig. ^ 
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Figure 6: Vt is plotted as a function of for the instances of Fig. ^ 
4.2 The diluted bimodal distribution 

This is the distribution in which the random variable defined by ()24() is replaced, for some 
values chosen at random, by 0. The resulting random variable is 

z = uy (26) 

where y is distributed according to 1)24(1 and u is distributed according to 

Piu) = p5il - u) + (1 - p)6{u) , (27) 

where < p < 1 is the dilution parameter. Thus, in our notations, p = 1 corresponds to 
no dilution. The mean of the distribution of z is 0, and its variance is p. We again apply 
the regularization process to the matrix A by multiplying the matrix elements Aij, chosen 
at random in this ensemble by ()25|) . 

The calculations were repeated for p = 0.5 and p = 0.2. In all cases a scaling function 
was found. The distributions Pa as a function of x'^ are presented in Fig. Qfor p = 0.5, 
and in Fig. Elfor p = 0.2. The analogous graphs of the distribution functions Pjs and Pt 
for these diluted bimodal ensembles are presented in |^. Figs. |7|and|Sl as well as the 
corresponding figures for Pp and Pt indicate that the convergence in the limit m — > oo 
becomes slower as the dilution increases (i.e., as p decreases). 
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Figure 7: Pa is plotted as a function of x'^ for the diluted bimodal distribution, where 
p = 0.5. As before, n = 2m. The number of instances used in the simulation is 54951 for 
the m = 20 case, 41107 for the m = 30 case and 50863 for the m = 40 case. The number 
of converging instances for each case is 5000. 
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Figure 8: Pa is plotted as a function of x'^ for the diluted bimodal distribution, where 
p = 0.2. As before, n = 2m. The number of instances used in the simulation is 54620 for 
the m = 20 case, 37697 for the m = 30 case, and 65367 for the m = 40 case. The number 
of converging instances for these cases were, respectively, 4980, 3725 and 5921. 
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Figure 9: Va is plotted as a function of x'^ for the uniform distribution. As before, n = 2m. 
The number of instances used in the simulation is 121939 for the m = 20 case, 91977 for 
the m = 30 case and 112206 for the m = 40 case. The number of converging instances for 
each case is 20000. 



4.3 The uniform distribution 

The distribution of the elements of j4, 6 and c in this ensemble is 

1 -l/2<y<l/2 

otherwise 



Piy) 



(28) 



For this distribution the mean is and the variance is j^- The distribution function H18() is 
presented in terms of the scaling variable (|2H) in Fig. El All the other distributions related 
to this ensemble are presented elsewhere j26l| . 
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Figure 10: Pa as a function of xa (for n = 2m and m = 40). The graphs are scaled to fit 
the theoretical Gaussian result by appropriate choice of the factors a!~^ . 

5 Universality 

In the previous section we demonstrated that for large m and n, the distribution functions 
(fT8|) - (|20|) depend on A,I3,T and on m and n only via the scaling variables (|2T|) - (|23|) . 
A natural question which arises is whether the distribution functions Paj P/3 and Pt are 
universal |12[ ll3j. In other words, we ask whether all probability ensembles of LP problems, 
or at least a large family thereof, yield the same functions Paj P/3 and Pt of the scaling 
variables 

XA = a^^\n/m)x'^, (29) 
xp = a'p\n/m)xp (30) 

and 

XT = a!^\n/m)x'rp (31) 

where x'^, x'^ and x'j, are defined by l|21|l . (|22|) and (|23j) . and the scale factors are the 

generalizations of a^^^ and Oj) of (fT^ . ((T3|) and (fTB|) ? 

For the Gaussian distribution we denote jj, = g, for the undiluted bimodal distribution 
/X = p = 1, for the various diluted bimodal distributions n = p := 0.5,0.2, and for the 
uniform distribution fi = u. Throughout this paper we take n = 2m. The distributions of 
the 

^min for m — 40 for the various ensembles are presented in Fig. 1101 
The scale factors are chosen so as to minimize the deviation of the specific distri- 
bution Pa from the Gaussian distribution J^{xa ) of Q , that was calculated analytically in 
^lE]. This is done, as usual, by least squares fit. The Gaussian distribution was taken 
with variance af^ = 1 in our calculation. For the Gaussian distribution the numerical 
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results were found to fit the analytical result © with ~ 0.564 ~ as expected from 
H12|) for n = 2m. For the bimodal distribution it was found that a^^ ~ 0.558 ~ For 
the diluted bimodal distributions at p = 0.5 we found a^'^^ ~ 0.581 ~ "^'^ ^ for the 
more diluted ensemble at p = 0.2 we found a^'^^ ~ 0.687 ~ ^ • Finally, for the uniform 

distribution we obtained a"^ ~ 1.957 ~ ~ ^J^^^- . Recall that the variances of these 

distributions are cr^^^ = 1, and ct^^^ = ^j^, respectively. Therefore, on the basis of these 
numerical results we conjecture that 

id) (1) (") /ori\ 

^(9)«A = ^(i)«A = '^{u)a\ ' (32) 

and that their common value is given by ()12() . that is equal to l/-v/vr for n = 2m. 

(p) 

Our numerical results indicate that the scale factors a\ of the diluted distributions 
deviate from this simple law, and this deviation seems to be more pronounced for higher 
dilution (smaller p). For these distributions cr^p^ = p, leading to ^^^^(o.s) — 0.411 and 

a^'^"* £7(0.2) — 0.307, which differ significantly from the more or less common value of this 
product for the undiluted distributions, namely 1/\/tt ^ 0.564. 

From Fig. ^Jwe see that the distribution functions of the scaling variables corre- 
sponding to the convergence rates A approach a universal function, that is identical to the 
one that is found analytically for the Gaussian distribution, and is given by @. For the 
undiluted distributions also the scale factors were found to agree with (|12j) . 



The proportionality of to l/c(^) probably results from the fact that if all parameters 
Cj and Aij in (j2} are rescaled by some common factor, also the Aj, defined by ©, are 
rescaled by the same factor. For the diluted distributions, the behavior of the scaling 
factors is different. Behavior of similar nature is found also for the distribution functions 
P/3 and Pt. 

The distributions Pp in terms of are presented in Fig. ^2 for m = 40. Using the 
scale factor = 1, we find numerically, by least square fit to the case of the Gaussian 
distribution (that was also found numerically), the other scale factors to be a^p ~ 0.952, 
a^^-^^ ~ 1.189, af^^ ~ 1.943 and aj^"^ ~ 0.960. 

The scale factors a^^\a^p^ and a^"'* are pretty close to each other (and to unity), while 

a^^'^'^ and a^*^'^^ deviate from them significantly. 

We note from that the barriers Pi are invariant under global rescaling of all the 
matrix elements Aij by the same factor. Thus, as long as the optimal vertex x*^ does not 
have too many anomalously small or large components, Pp is expected to be independent 
of the variance, and the numerical values we obtained for a^p ~ a^"^ ~ a^^^ seem to support 
this expectation for the undiluted ensembles. The scale factors a^^'^\ a^^'"^^ deviate from 

that common value, with a^'^'^^ deviating very significantly. 

Our numerical results, displayed in Fig. show that the distributions Pp{xp) found for 
the diluted ensembles deviate from the ones found for the undiluted ensembles. Therefore, 
for the diluted ensembles, although we have good indication for universality, this point 
clearly requiers further research. 

The distributions Pt in terms of xt are presented in Fig. I12| for m = 40. 
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Figure 11: as a function of Xj^ for all the distributions checked (for n = 2m and 

m = 40), where the scale factors a^^^ were found by least squares fit to the distribution for 
the Gaussian ensemble, which was found numerically as well. 
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Figure 12: Pt as a function of for all the distributions checked (for n = 2m and m = 40). 
The scale factors Oj^^ were found by least squares fit to the distribution for the Gaussian 
ensemble, which was found numerically as well. 
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Using the scale factor aif^ = 1 we find numerically a^^ ~ 1.008, a^'^^ ~ 1.240, a^'^"* ~ 

2.088, and aj! '■' ~ 3.399 ~ Vll-553. From © we see that under global rescaling of the 
parameters Aij,bi,Ci in ©, the scaled variables x'^ and x'rp are rescaled in the same way. 

Thus, since is proportional to 1/cr, also a^'* should satisfy a similar proportionality. 

Since for the Gaussian distribution we have taken cr^^-j = 1 and = 1 , our numerical 

results, where we find ~ \/T2 , suggest that for the undiluted ensembles 

(s) (1) (") /oo\ 

^ia)"'T = f^(i)«T = ^^(m)«t ' (33) 

taking the value unity in our case. For the diluted ensembles we find o"(o.5)Q^t^^ — 0.877 and 

cr(0.2)Oj?^^ — 0.933, deviating significantly from unity. For p = 0.2 a significant deviation 
of Pt from the other distributions is found. 



6 Summary and Discussion 

In this paper we have presented ample numerical evidence for the fact that the asymptotic 
distribution functions Pa,P/3 and Pt are scaling functions. 

In particular, all the undiluted ensembles of LP problems which we studied, seem to 
have the same set of asymptotic distribution functions, when the latter are expressed in 
terms of the scaling variables XA,a^/3 and xt- Furthermore, we have found that then the 
following combinations of scaling factors, o"(^)a^'*, a^^^ and fT(^)a^^ are independent of the 
probability distribution. Therefore, in particular, should satisfy H12|l that was found 
for the Gaussian distribution, while a^^^^ and well as corresponding distributions, 

are yet to be found analytically. 

Based the results presented in this paper, as well as the results of ^[E], we conjecture 
the scaling behavior of the various undiluted distribution functions (and the corresponding 
scaling factors) is universal, i.e., that it is robust and should be valid in a large class of 
ensembles of LP problems, in which the {A, b, c) data are taken from a distribution with 
zero mean and finite variance. 

We have also studied the effect of dilution, namely, imposing that the parameters Aij,Ci 
and hi in (j2} could vanish with finite probability 1 — p. Specifically, we have studied the 
effects of dilution only on the bimodal distribution. 

Our findings, depicted in Figs. [3 and |H1 as well as the corresponding figures for 
and Pt, indicate that the convergence in the limit m — > oo becomes slower as the dilution 
increases (i.e., as p decreases). More importantly, our numerical results indicate that the 
diluted distributions may exhibit scaling behavior as well, but with scaling factors which are 
different from those of the undiluted ensembles which belong in the universality class of the 
undiluted Gaussian ensemble. Moreover, the corresponding asymptotic scaling distribution 
functions Pa , P/3 and Pt of these diluted ensembles deviate sometimes from those of the 
corresponding ones of the undiluted Gaussian universality class. This deviation appears 
to become more pronounced as dilution increases, and it may be related to the fact that 
a finite fraction of the admissible LP instances in diluted ensemble may have an optimal 
vertex x*^ which has too many anomalously small or large components. 
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Generally speaking, dilution seems to have interesting effects, which are not completely 
understood, and call for further investigation. Specific questions arc motivated by the 
present work. In particular, arc the asymptotic scaling distribution functions P^,Pp and 
Pt, which we computed numerically, really different from the ones found for the undiluted 
Gaussian ensemble (with possibly scale factors which differ from the Gaussian ones), or 
the deviations of these functions, indicated by our numerical results, are merely effects of 
the slower convergence towards asymptotics? If they are different - do they form another 
universality class? 

Finally, we would like to raise a question which may be of practical importance. Thus, 
imagine that all LP problems used in practice (or at least, a large fraction thereof) are 
collected into an unbiased probability ensemble. How is this distribution of realistic LP 
problems related to the ensembles studied in this paper? Does it really relate to a univer- 
sality class (or classes) of ensembles of LP problems studied here? Does it agree more with 
the diluted or the undiluted ensembles? These questions clearly pose important conceptual 
challenges for further investigation, and also have practical implications. 
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